clear;
clc;

%% Setup Information Structure
Info            =   struct();
Info.ParamList  =   {'KAPPA','b_lb','B_Total_SS'};
Info.ParamNum   =   length(Info.ParamList);
Info.ValueList  =   struct('KAPPA',(0:0.01:0.20)'/4,...
                           'b_lb',(-4:0.2:0)',...
                           'B_Total_SS',(0.1:0.2:4)');
Info.Value      =   struct('KAPPA',0.08'/4,'b_lb',-0.3,'B_Total_SS',1);
Info.Num        =   0;
UnitValue       =   zeros(1,Info.ParamNum);
for ii=1:Info.ParamNum
    vv              =   Info.ParamList{ii};
    Info.Idx.(vv)   =   Info.Num+(1:1:length(Info.ValueList.(vv)))';
    Info.Num        =   Info.Num+length(Info.ValueList.(vv));
    UnitValue(ii)   =   Info.Value.(vv);
end
Info.UnitValue  =   UnitValue;
ParamMat        =   repmat(UnitValue,[Info.Num,1]);

for ii=1:Info.ParamNum
    vv              =   Info.ParamList{ii};
    ParamMat(Info.Idx.(vv),ii) ...
                    =   Info.ValueList.(vv);
end

Info.Node       =   ParamMat;
Info.Num        =   size(Info.Node,1);
Info.StatName   =   {'Prob. of Negative Wealth',...
                     'Median Wealth/Median Income',...
                     'Wealth-Income Ratio: Median',...
                     'Prob. of HtM',...
                     'MPC: Median','MPC: Average' ...
                     };
%% Solve the Moments at Different Parameter Values
ResultCell      =   cell(Info.Num,1);
TempFun         =   @(ParamVec)Study_1_SubFun_TargetMom(Info,ParamVec);
parfor ii=1:Info.Num
    ResultCell{ii}  =   Study_1_SubFun_TargetMom(Info,ParamMat(ii,:));
end
BenchmarkResult     =   Study_1_SubFun_TargetMom(Info,UnitValue);
save('Identification.mat','Info','ResultCell','BenchmarkResult');

%% Checks
TempStatIdx     =   [4;5];
NumStat         =   length(TempStatIdx);
Moment          =   zeros(NumStat,Info.Num);
parfor ii=1:Info.Num
    if ResultCell{ii}.Flag
        Moment(:,ii)    =   ResultCell{ii}.MomVec(TempStatIdx);
    else
        Moment(:,ii)    =   NaN(NumStat,1);
    end
end

BenchmarkMoment     =   BenchmarkResult.MomVec(TempStatIdx);
%% Figure
Folder      =   'TableGraphs/Identification/';
mkdir(Folder);
ParamDict   =   struct('b_lb','$\underline{b}$','B_Total_SS','$\bar{B}$');
ParamList   =   fieldnames(ParamDict);
StatList    =   {'Median Wealth-Income Ratio','Median MPC'};

N_col       =   length(ParamList);
N_row       =   length(TempStatIdx);
fig_size    =   [N_col,N_row].*[1,1]*1/4; 
fig_gap     =   [0.10,0.10];
fig_VMargin =   [0.05,0.05];
fig_HMargin =   [0.10,0.05]; 

fig_SubPlot =   @(ii_vv)subtightplot(N_row,N_col,ii_vv,...
                                     fig_gap,fig_VMargin,fig_HMargin);

Fig         =   FigureSetup(['Steady State Moments'],fig_size);
for jj=1:N_row
    for kk=1:N_col
         ax     =   fig_SubPlot((jj-1)*N_col+kk);
         vv     =   ParamList{kk};
         plot(Info.ValueList.(vv)',Moment(jj,Info.Idx.(vv)),'-b','LineWidth',2);
         hold on;
         yline(BenchmarkMoment(jj),':k','LineWidth',1);
         xline(Info.Value.(vv),':k','LineWidth',1);
         ylim([0,0.8]);
         title(ParamDict.(vv),'Interpreter','Latex','Fontsize',8);
         ylabel(StatList{jj},'Interpreter','Latex','Fontsize',8);
         ax.YGrid   =   'on';
         ax.XGrid   =   'on';
    end
end
print('-depsc','-r1000',Fig,[Folder,'WealthIncomeDistribution']);